Chaotic orbits in thermal-equilibrium beams: existence and 

dynamical implications 

Courtlandt L. Bohn 1,2 and Ioannis V. Sideris 1 
1 Northern Illinois University, 

DeKalb, IL 60115 
2 'Fermilab , Batavia, IL 60115 
(Dated: February 2, 2008) 

Abstract 

Phase mixing of chaotic orbits exponentially distributes these orbits through their accessible 
phase space. This phenomenon, commonly called "chaotic mixing", stands in marked contrast 
to phase mixing of regular orbits which proceeds as a power law in time. It is operationally ir- 
reversible; hence, its associated e-folding time scale sets a condition on any process envisioned 
for emittance compensation. A key question is whether beams can support chaotic orbits, and if 
so, under what conditions? We numerically investigate the parameter space of three-dimensional 
thermal-equilibrium beams with space charge, confined by linear external focusing forces, to de- 
termine whether the associated potentials support chaotic orbits. We find that a large subset of 
the parameter space does support chaos and, in turn, chaotic mixing. Details and implications are 
enumerated. 

PACS numbers: 41.75.-i, 05.70.Ln, 29.27.Bd, 45.10.Na, 98.10,+z 
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I. INTRODUCTION 



Rapid, inherently irreversible dynamics is a practical concern in producing high-brightness 
charged-particle beams. Time scales of irreversible processes place constraints on methods 
for compensating against degradation of beam quality caused by, for example, space charge. 
This is a very important practical matter because compensation must be fast compared to 
these processes, and this affects the choice and configuration of the associated hardware. 

A beam bunch with space charge comprises an iV-body system with typically 3N degrees 
of freedom. Upon coarse-graining, i.e., "smoothing" the system to remove granularity, the 
collective space-charge force remains. One might conjecture that this force, when nonlinear, 
may support chaotic orbits. One example is the University of Maryland five-beamlet experi- 
ment that shows presumably irreversible dissipation of the beamlets after a few space-charge- 
depressed betatron periods Simulations of the experiment reveal a substantial fraction 
of globally chaotic orbits Q], and phase mixing of these orbits thereby presents itself as a 
contributing evolutionary mechanism. This example pertains to a strongly time-dependent 
nonequilibrium system, yet one might conjecture that nonlinear space-charge forces in a 
static system could support chaotic orbits as well. We shall explore this conjecture. 

An initially localized clump of chaotic orbits will, via phase mixing, grow exponentially 
and eventually reach an invariant distribution. This is "chaotic mixing" ^, Q]. Strictly 
speaking, the process is reversible in that it is collisionless and its dynamics is included in, 
e.g., Vlasov's equation. Nonetheless, when the invariant distribution spans a global region 
of the system's phase space, chaotic mixing is a legitimate relaxation mechanism in that it 
drastically smears correlations. Moreover, from a practical perspective, the process is strictly 
irreversible because infinitesimal fine-tuning is needed to reassemble the initial conditions. It 
is also distinctly different from phase mixing of regular orbits, i.e., linear Landau damping j^|, 
a process that winds an initially localized clump into a filament over a comparatively narrow 
region of phase space. Whereas chaotic mixing proceeds exponentially over a well-defined 
time scale and can cause global, macroscopic changes in the system, phase mixing of regular 
orbits carries a power-law time dependence, proceeds on a time scale depending on the 
distribution of orbital frequencies across the clump, and acts only over a portion of the 
phase space. Accordingly, ascertaining conditions for, and time scales of, chaotic mixing in 
beams is an undertaking of practical importance. 
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In this paper we consider a family of thermal-equilibrium (TE) configurations of 



with the bunch and has its origin affixed to the bunch centroid. Particle motion in this 
reference frame is taken to be nonrelativistic; transforming from the bunch frame to the 
laboratory frame is straightforward [9J. In the laboratory frame the space-charge force de- 
creases inversely with the square of the beam energy. For the transverse component, this 
arises from the partial cancellation between the self-magnetic and self-electrostatic forces; 

rn 

while for the longitudinal component, it is due to Lorentz contraction jlOl ]. Nonetheless, 
there are many situations involving high-brightness beams wherein space charge is impor- 
tant. Contemporary examples include low-to-medium-energy hadron accelerators such as 
those that drive spallation-neutron sources or serve as boosters for high-energy machines, 
heavy- ion accelerators, and low-energy electron accelerators such as photoinjectors 



Thermal-equilibrium beams are of practical interest in connection with, e.g., high-current 
radiofrequency linear accelerators. While conventional designs of such machines lead to 
bunches that are out of equilibrium, a design strategy that keeps the beam at or near 
thermal equilibrium has been formulated [12]. The principal motivation for this alternative 
strategy is to circumvent equipartitioning processes that cause emittance growth and halo 
formation. 

Because a TE configuration is a maximum-entropy configuration, is static, and is mani- 
festly stable , one might expect its intrinsic dynamics to be entirely benign. The expecta- 
tion is questionable. The density distribution of such a configuration is uniform in its interior 
and falls to zero over a distance commensurate to the Debye length. Thus, large-amplitude 
orbits will explore this "Debye tail," during which time they experience a nonlinear force. 
The question we seek to answer is whether the nonlinear force in the Debye tail can cause a 
significant number of orbits to be chaotic. The answer is unequivocally "no" for spherically 
symmetric or infinitely long cylindrically symmetric configurations because their potentials 
are integrable and thereby support only regular orbits. However, breaking the symmetry 
can generate chaotic orbits, as will become apparent in the analysis to follow. 

Our study involves a comprehensive suite of numerical experiments concerning orbital 
dynamics in smooth (coarse-grained) TE configurations. We establish a quantitative mea- 
sure of chaos in orbits and use this measure to distinguish between regular and chaotic 
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orbits. We then evolve initially localized clumps of particles in the smooth potentials. The 
experiments are fast if the potentials are analytic, but they are much slower if the potentials 
must first be tabulated numerically over a grid. As part of the preliminaries, Sec. HTlpresents 
a semianalytic theory for estimating the time scale for chaotic mixing. In general the TE 
configurations, specified in Sec. Illll must be found numerically. Section |TV| presents a means 
for rapidly constructing approximate, semianalytic models of their potentials. With these 
models we are able to survey the parameter space and obtain a zeroth-order assessment of 
the prevalence and degree of chaos; this is done in Sec. |V] Section IVT1 concerns examples for 
which the potential is accurately determined via a numerical solution of Poisson's equation 
on a grid. For these examples the experiments of Sec. |3 are repeated, and the results are 
compared to those derived from the semianalytic approximation. Section IV 111 summarizes 
the findings, discusses their implications while providing a comparison with the theory of 
Sec. m and presents a path for follow-on work. 



II. ESTIMATED TIME SCALE FOR CHAOTIC MIXING 



Before embarking on numerical studies, it is wise to ascertain whether chaotic mixing can 
indeed proceed rapidly. One can construct an analytic tool to estimate the chaotic-mixing 
rate, although its application involves the tacit assumption, or initial knowledge, that chaotic 
orbits are present. In this section we sketch the methodology leading to analytic predictions. 



Additional details are available elsewhere 



14. 



The past few years have seen development of a geometric method proposed by M. Pettini 
to quantify chaotic instability in Hamiltonian systems with many degrees of freedom. The 
central idea is to describe the dynamics in terms of average curvature properties of the 
manifold in which the particle orbits are geodesies. The method hinges on the following 
assumptions and approximations; they are discussed thoroughly in Ref. (1) a generic 
geodesic is chaotic; (2) the manifold's effective curvature is locally deformed but otherwise 
constant; (3) the effective curvature reflects a gaussian stochastic process; and (4) long- 
time-averaged properties of the curvature are calculable as phase-space averages over an 
invariant measure, specifically, the microcanonical ensemble. The gaussian process is the 
zeroth-order term in a cumulant expansion of the actual stochastic process; assumption 
(3) is that the zeroth-order term suffices. The end result relates chaotic instability to the 
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geometric properties of the manifold defined by the long-time-averaged orbits. In short, the 
theory is based on (often questionable) assumptions that chaos exists and is characterized 
by ergodicity and a microcanonical ensemble, and it treats chaotic orbits as arising from 
a parametric instability that can be modeled by a stochastic-oscillator equation. It has 
recently been adapted for application to low- dimensional, autonomous (time-independent) 
Hamiltonian systems and, in tests against a wide variety of such systems, it was found 
commonly to yield estimates of mixing rates that are good to within a factor ~ 2 jis| . 

Action principles in classical mechanics are tantamount to extremals of "arc lengths"; 
thus, one can infer a m e,ie tensor fr0m an action princip.0 Q. The metric ten 80r ,nan lfe8 t 8 
all of the properties of the manifold over which the system evolves, with these properties 
being calculable following standard methods of differential geometry. Of special interest is 
the divergence of two initially nearby 3iV-dimensional geodesies q and q + 5q as governed 
by the equation of geodesic deviation: 

-^-+R Pjt —6q — - t (1) 

in which D/ds denotes covariant differentiation with respect to the "proper time" s, R a /3js is 
the Riemann tensor derivable from the metric tensor, and summation over repeated indices 
is implied with each index spanning the 3N degrees of freedom. Equation is the basis 
for determining the mixing rate \ as a measure of the system's largest Lyapunov exponent, 
a quantity that reflects the long-time behavior of the separation vector: 

x = IhnilnMt. (2) 



Any number of action principles, and therefore any number of metric tensors, can 
be selected to proceed further. Eisenhart's metric [18], which is consistent with Hamil- 
ton's least-action principle, is probably the most convenient choice. It offers easy cal- 
culation of the Riemann tensor, and it avoids spurious results traceable to the singu- 
lar boundary of the perhaps better-known Jacobi metric that is derivable from Mauper- 
tius' least-action principle [3]. Eisenhart's metric operates over an enlarged configura- 
tion space-time manifold in which the geodesies are parameterized by the real time t, i.e., 
ds 2 = dt 2 = — 2V(q)(dq ) 2 + bijdtfdqi + 2dq°dq 3N+1 , in which V(q) is the potential energy 
per unit mass (hereafter called the "potential"); 5ij (with the indices i, j running from 1 to 
3N) is the unit tensor corresponding (without loss of generality) to a cartesian spatial coor- 
dinate system, q° = t; q 3N+1 = t/2 — f c?t'L(q, q); and L is the Lagrangian. The resulting 



geodesic equations for the spatial coordinates q l are Newton's equations of motion, so the 
particle trajectories correspond to a canonical projection of the Eisenhart geodesies onto the 
configuration space-time manifold. A convenient byproduct of the Eisenhart metric is that 
the only nonzero components of the Riemann tensor are -RoiOj = didjV, in which <9j — 91 dq l . 

Using the aforementioned assumptions and approximations, Pettini and others jl6L |20( 
derive an expression for \ in terms of the curvature and its standard deviation averaged over 
the microcanonical ensemble. The idea is that, as t — > oo, chaotic orbits of total energy E 
mix through the configuration space toward an invariant measure, taken per assumption (4) 
to be the microcanonical ensemble 5(H — E), over which time averages become equivalent 
to phase-space averages. Specifically, for an arbitrary function A(q), the averaging process 
is 

Per Eisenhart's metric, the average curvature k and the ratio p = a/n, with a denoting the 
standard deviation of the curvature, are 



(v 2 v) _ i vav 2 ^) 2 -^) 2 ) f . 

K 3N-V 9 k vW^I ' {) 

in which V* denotes the Laplacian 00, and p corresponds physically to the ratio of the 
average curvature radius to the length scale of fluctuations [21]. By taking the curvature 
to vary randomly along a chaotic orbit, one can reduce Eq. ((TJ to a stochastic-oscillator 
equation that can be solved analytically. The solution yields an estimate of the largest 
Lyapunov exponent \'- 

1 L\p) - 1 



x(p) 



L(P) 



VS L(p) 



T(p) + Jl + T*(p) 



' T{P) = 8 2VTT-p + np (5) 



The geometric quantities derive from the 6A^-dimensional microcanonical distribution. 
Anticipating that granularity takes a long time to affect mixing, and wishing to identify con- 
ditions for rapid mixing, we now consider the influence of the 3-dimensional coarse-grained 
space-charge potential V s on a generic chaotic orbit. The largest Lyapunov exponent for the 
coarse-grained system equates to the chaotic-mixing rate. We presume the assumptions and 
approximations stated at the outset carry over to the coarse-grained system; the main jus- 
tification is that the aforementioned previous work concerning low- dimensional autonomous 



Hamiltonians has shown the mixing rate in such systems usually depends only weakly on 
the dynamical details We take the external focusing potential Vf to be quadratic in the 
coordinates x comoving with the bunch, i.e., V/(x) = (u ■ x) 2 /2, wherein u = (aj x ,(j y ,u z ) 
corresponds to the focusing strength; the total potential is V — Vf + V s . Per Eq. (0J) and 



are the single-particle charge and rest mass, respectively, and e Q is the permittivity of free 
space. We then have 



Inserting these results into Eq. @ gives the associated time scale for chaotic mixing, 
t m = 1/x- When the standard deviation of the density distribution is large, as can be 
the case when substructure is present, p will be appreciable, and in turn Eq. (JSJ) makes 
clear that t m will be a few space-charge-depressed periods It^I^Jk. Accordingly, the space- 
charge-depressed period, a quantity commensurate to the orbital period of a typical particle, 
constitutes a "dynamical time" tu for charged-particle beams. 

To underscore the potential impact of collision /ess relaxation via chaotic mixing, it is 
of interest to compare t m to the collision al relaxation time t R . Perhaps the simplest way 
to develop an order-of-magnitude estimate of t R in a charged-particle bunch (a nonneutral 
plasma) is to calculate the time required for a typical particle velocity to ch ange by of order 
itself presuming collisions comprise a sum of incoherent binary interactions [22]. The result 
is tn/tp ~ O.liV/lniV, wherein the Coulomb logarithm is conservatively taken to be IniV. If 
we substitute plausible parameter values for real high-brightness beams, we find t R ^> to; 
for example, N = 6.25 x 10 9 (1 nC) gives t R ~ 10 7 tp; hence, t m <C t R when chaotic mixing 
is prominent. The remaining question is whether there is a significant population of globally 
chaotic orbits to mix, a question to which we now turn our attention. 

III. THE EQUATIONS OF THERMAL EQUILIBRIUM 

Consider a system, i.e., a bunch, of N identical charged particles, e.g., electrons or pro- 
tons. For simplicity, invoke a Cartesian coordinate system whose origin lies at the bunch 
centroid. Assume all particle velocities in this coordinate system are nonrelativistic. The 



Poisson's equation the quantities k and a are determined from V 2 V = u"j — Wp(x), in which 
uPi = ul + Uy + oj 2 z1 ^p(x) = n(x)g 2 /(e m), n(x) is the (smoothed) particle density, q and m 
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particles mutually interact via the Coulomb force and are confined by a static, externally 
applied, linear focusing force. The focusing force may have different strengths along the 
three Cartesian axes. Assume, apart from this focusing force, that the system is isolated 
and is in thermal equilibrium. Accordingly, the total energy E of each particle is conserved: 

E = \-rnv 2 + ]rm(uj ■ x) 2 + q<j)(x) ; (7) 

wherein to = (a> x , u y , u z ) corresponds to the focusing strength; x = (x,y,z) denotes co- 
ordinates; m, v, and q are the particle's rest mass, speed, and charge, respectively; and 
0(x) = (m/q)V s is the space-charge potential arising from the collective Coulomb force. 

To proceed, one would in principle work with the 6iV-dimensional microcanonical dis- 
tribution of particles. This distribution includes interactions at all scales, ranging from 
particle-on-particle to a single particle interacting with the bulk, smooth potential from all 
other particles. Discreteness effects from 1/r 2 particle collisions generate chaos j^J; they 
cause nearby particle trajectories to separate exponentially. The rate of exponential sep- 



aration, i.e., the Lyapunov exponent, is an increasing function of iV 24]. In this sense, 
larger N gives rise to more chaos. However, the scale at which the separation saturates is 
a decreasing function of N. Accordingly, in large-iV, high- charge- density systems such as 



beams with space charge, discreteness establishes microchaos [25|, |26|, |27|, |28J . At the other 
extreme, that of a single particle interacting with the bulk, smooth potential, exponential 
separation of nearby chaotic particles (if any are present) saturates at a global scale, cor- 
responding to a state of macrochaos. Thus, initially nearby chaotic orbits evolve in three 
stages 2^1 : (1) very rapid exponential divergence that saturates at a scale large compared to 
the initial interparticle spacing but small compared to the system size; followed by (2) rapid 
exponential divergence that persists until the particles are globally dispersed; followed by 
(3) less rapid power-law divergence on a time scale oc (In N)to, in which to is a dynamical 
time commensurate to the orbital period. If, in the smooth potential, the initially nearby 
particles execute regular motion rather than chaotic, then stage (2) is absent, and stage (3) 
proceeds on the much longer time scale cx (N 1 / 2 )tr> 27 1. 

Our interest here is in stage (2). Specifically, we are concerned about the existence of, and 
time scale for, macroscopic chaos, i.e., chaotic mixing into the global region of phase space 
that is energetically accessible to the individual particles. Accordingly, we specialize to the 
smooth 6-dimensional distribution function of a single particle, recognizing that discreteness 
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effects vanish on macroscopic scales as the number density grows. For the TE beam, this 
is just the Maxwell-Boltzmann distribution, /(x, v) oc exp(—H/kT), in which H = E is 
the Hamiltonian, k is Boltzmann's constant, and T is the beam temperature. The number 
density follows upon integrating over velocity space, and the space-charge potential follows 
upon solving Poisson's equation: 



n(x) = n(0) exp 



\m(oj ■ x) 2 — g0(x) 



kT 



(8) 



V 2 0(x) = --n(x), 0(x = 0) = V0(x = 0) = 0, (9) 

wherein e Q is the permittivity of free space. 

A much more convenient formulation arises by using dimensionless variables. We intro- 
duce the Debye length X D0 and angular plasma frequency u p0 , both defined in terms of the 
centroid number density n(0): 

2 _ e kT 2 _ n(0)q 2 
n(U)q z F e m 

We then measure all lengths in the unit of A^o, i.e., x «-> x/A/xb and all times in the 
unit of 1/ujpo, i.e., t <-> o; p o£- I n addition, we introduce the dimensionless potential <&(x) = 
q4>(x) / (kT) , and we normalize n(x) to the centroid density n(0), i.e., n(x) <-> n(x)/n(0). 
The number density and Poisson's equation then reduce to their dimensionless forms: 

1 



n{x) = exp 



2 -fi 2 i? 2 (x) - $(x) 



(11) 



V 2 $(x) = -n(x), $(x = 0) = V$(x = 0) = 0; (12) 

wherein fi 2 = (cjy/cjpo) 2 ^ and -R 2 (x) = (x/a) 2 + y 2 + (z/c) 2 , with the "scale lengths" 
a and c defined as a = uj y /uj x and c = io y /u z . 

Equations and self-consistently provide the structure of the entire family of 
smoothed TE configurations. The parameter Q governs the strength of external focusing vis- 
a-vis space charge. The scale lengths a and c set the overall geometry: a = c = 1 corresponds 
to spherical symmetry, a = c ^ 1 corresponds to cylindrical symmetry, and a ^ 1, c ^ 1, 
a 7^ c establishes a triaxial configuration. The extreme case of maximum space charge 
corresponds to a density that is strictly uniform over the volume of the configuration, in 
which case $ = —Q 2 R 2 /2 inside the beam. Upon substituting into Poisson's equation, we see 



that the associated configuration carries the parameter = fi u = 1/ J(l/a 2 ) + 1 + (1/c 2 ). 



This is the minimum permissible focusing strength; the bunch is unconfmed if Q < Q u , and 
the corresponding constraint on the parameter space is 

1 1 1-Q 2 

< 13 > 

Hence, the parameter set [a, c; Q] fully specifies a TE configuration. 

Upon solving for the space-charge potential $(x), one can calculate orbits of test particles 
in the total potential. Their trajectories follow from the (dimensionless) equation of motion: 



d 2 x 



(Jf2 -V [-fi 2 i? 2 (x) + $(x)j . (14) 

One can, of course, introduce arbitrary initial conditions for the orbits. In our experiments, 
the initial condition on the velocity is v(0) = 0, and the total energy E of a particle thereby 
corresponds to the potential energy associated with the initial position x(0). 

A key challenge in exploring orbital dynamics throughout the parameter space is to 
integrate large numbers of orbits rapidly for sufficiently long evolutionary times. Ideally, 
one would have analytic solutions for the density-potential pairs, from which the force on 
a particle at each time step can be quickly evaluated. Unfortunately, the equations of 
equilibrium generally do not submit to analytic techniques. Thus, in principle, one must 
solve these equations numerically, e.g., over a grid. However, as delineated in the following 
section, it is possible to formulate approximate, semianalytic solutions, and these solutions 
enable a search of a broad range of the parameter space for regions that support chaotic 
orbits. We now turn to that exploration. Subsequently, for select cases, we compare these 
results against those derived from fully self-consistent numerical solutions. 



IV. APPROXIMATE SOLUTIONS TO THE EQUATIONS OF EQUILIBRIUM 

A method to solve the equations of equilibrium is through a sequence of successive ap- 
proximations [3^. A way to begin such a sequence is as follows: (1) As a first approximation, 
represent the system as a configuration stratified on similar and similarly situated concen- 
tric ellipsoids. A "homeoid" is a shell that is bounded by two similar and similarly situated 
concentric ellipsoids, and in which the surfaces of constant density are ellipsoids that are 
similar to, concentric with, and similarly situated with respect to the bounding ellipsoid. 
Thus, the charge density is "homeoidally striated" in the first approximation (as is later il- 
lustrated in Fig. IT^|). Determine the stratification by solving a spherically symmetric model 
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of the equations of equilibrium. (2) In the second approximation, derive the space-charge 
field corresponding to the homeoidally striated charge density, and then solve exactly the 
equations of equilibrium in this field. (3 and up) Repeat the process until the density and 
potential converge. In practice, one can carry out steps (1) and (2) of this recipe using 
semianalytic methods; to go further requires numerical techniques. 



A. Determination of the structure in the first approximation 



To invoke a spherically symmetric model of Eq. (|12|). we take the potential to be stratified 
over ellipsoids on which i?(x) takes a constant value. Then the spherically symmetric model 
corresponds to solving 

1 



1 d 
R^dR 



$o(0) 



— exp 
d$ 



dR 



n 2 R 2 - $ (R) 



0. 



(15) 



R=0 



This model defines the "zeroth approximation" <&o(R) to the potential. In general Eq. (JT5J) 
must be solved numerically; however, the solution is rapidly and easily accomplished with 
the aid of, e.g., a Runge-Kutta algorithm. 

Once &o(R) is determined, the corresponding homeoidally striated density becomes the 
first approximation to the number density: 

1 



ni(R) = exp 



~q 2 r 2 - $ (i?) 



(16) 



By inspection 



one can write down the space-charge potential corresponding to the 
number density ni(R), and this becomes the first approximation to the potential: 

ac r°° du r R ( x ; u ) 



$i(x) = "7T / XTT / drr ni (r) 
2 Jo A[u) jo 



ac 



du 



2 Jo A(u) 



dr 



(17) 



- r=R(x;u) 

wherein the second equality follows from an integration by parts, and the quantities A(u) 
and -R(x; u) are 



i?(x; u 



(u) = V '(a 2 + u){l + u){c 2 + u) 



a 2 + u 1 + u c 2 + u 
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Hence, in the first approximation the number density is homeoidally striated, but the space- 
charge potential is not. 

B. Determination of the structure in the second and higher approximations 

The number density in the second approximation, n2(x), follows upon substituting $i(x) 
calculated from Eq. (|17j) into Eq. (jllj) . For the special case of spherical symmetry, all orders 
of approximation agree with one another, but this is of course not true for a general triaxial 
geometry. To go further requires numerical methods, e.g., solving Poisson's equation for the 
potential $2 corresponding to the density n 2 , substituting the result into Eq. (fTTj) to obtain 
n 3 , and successively repeating the process until convergence is achieved. As discussed in 
Sec. IVII below, we use a different method, a multigrid algorithm, for solving Eqs. (jllj) and 
(|12|) numerically. 

V. SURVEY OF THE PARAMETER SPACE 

Gathering sufficient data to support precise, statistically based conclusions concerning 
orbital behavior in a given potential requires integrating thousands of orbits in that potential. 
And before these orbits can be tracked, the potential needs to be ascertained to sufficient 
accuracy. In principle, and for each choice of parameters, one must construct the "exact" 
potential $(x) by numerically solving the corresponding Poisson equation. This can be a 
computationally tedious process, and the solution is by necessity defined over a grid. Next, 
orbit integration through the grid requires accurate interpolation to evaluate the potential 
and corresponding particle acceleration between grid points. For sufficient resolution, the 
time steps need to be appropriately small; accordingly, many interpolations are required, 
and integrating many orbits is computationally time-consuming. This process is feasible 
for studying a few choices of parameter sets, and it underlies the results of Sec. IVII below. 
However, to survey the entire parameter space, i.e., to investigate many choices of parameter 
sets, the process becomes prohibitive. For this purpose one must resort to using approximate 
potentials. 

Sec. IIVI above details a sequence of approximations, the first elements of which are semi- 
analytic. The zeroth-order potential $o ; derived from Eq. (fTo^) . is easy to evaluate, and it 
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enables fast, high-precision orbital integration. However, $o itself may be a crude approxi- 
mation to the exact potential; the approximation gets progressively worse as the parameter 
sets deviate further from spherical symmetry. One might expect the potential $1 of the 
first approximation to provide a better model. However, its underlying integral, given in 
Eq. (|17|). adds additional complexity and time to the orbit integrations. We tried evaluating 
this integral at each time (thus position) step along the orbit, but doing so made the orbit 
computations prohibitively long. The alternative is to evaluate the integral over a grid and 
then do orbit integrations through the grid. As previously mentioned, integrations through 
a grid are too computationally expensive to enable a parameter survey. Moreover, if one 
is able to solve Poisson's equation for the exact potential $(x), then there is neither com- 
putational benefit nor motivation for using $1. Our strategy is to explore a few choices of 
parameter sets in the exact potential to strengthen conclusions from our survey, and this 
necessitated developing the Poisson solver described in Sec. IVII For these reasons, we use the 
potential of the zeroth approximation, <3>o, to survey the parameter space. For a few specific 
parameter sets for which the results of the zeroth approximation look especially interesting, 
we then check the results using the numerically evaluated exact potential, $(x). 



A. Solution for the zeroth-order potential 

Per Eq. ()1HJ) . the minimum possible focusing strength corresponding to a spherically 
symmetric system is £1 = Cl u = 1/ y/3. Accordingly, we choose focusing strengths in keeping 
with the following labeling convention: 



^7S 



1 + lO 1 ^ - , (18) 



with i = 1,2,..., I max and I max = 10. Choices of parameters [a, c; then must be selected 
based on the constraint of Eq. (|T3J) which bounds the parameter space of exact potentials. 
Note, however, that we can examine any desired geometry: oblate axisymmetric (for which 
a — 1, c < 1), prolate axisymmetric (for which a = c < 1), and the full range of oblate- 
through-prolate triaxial systems. Hereafter we refer to "Case 1, Case 2,..." according to 
"i = 1,2, ..." in Eq. (JUJ), respectively. 

Plots of <&o(R) versus R derived from Eq. (|15p appear in Fig. ^ Also shown are the 
corresponding profiles of the number densities n\(R) constructed in the first approximation. 



13 



For larger "case numbers" i, the density contains larger quasi-uniform central regions. In the 
outer regions the density decreases, over a length commensurate to the Debye length, to a 
low-density tail. The space-charge force in the quasi-uniform "core" is correspondingly quasi- 
linear; however, it is manifestly nonlinear in the "Debye fall-off region" (henceforth called 
the "Debye tail"). Fig.[T]shows that the choices of per Eq. (jTHj) span a wide range of space 
charge. At the one extreme, zero space charge, the density profile is gaussian. Then, the 
range of i spans from small space charge (i — 1) for which the density profile is approximately 
gaussian, through the fully space-charge-dominated, uniform beam {i = I max = 10) for which 
= l/y/3. Note that Case 5 (fi ~ l.OOOl/v^) represents "intermediate space charge"; the 
density falls off over a length scale comparable to that of the core. 



B. Methodology for orbital analysis 

1. Samples of orbits, power spectra, and complexity 

After choosing a parameter set [a, c; flj and solving for $ [-R(x)], we began by generating 
2000 initial coordinates uniformly spanning the volume occupied by the core and Debye tail 
of the density ni[R(x)]. The initial velocities were all chosen to be zero. Next, starting with 
these initial conditions, for every orbit we integrated the equations of motion, cf. Eq. (|14j) . 
using $ = $o for at least 100 orbital periods and, in most cases, for ^ 200 orbital periods. 
The integrations were done using a fifth-order Runge-Kutta algorithm with variable time 
step. As the integration proceeded, we computed the largest short-time Lyapunov exponent 
of each orbit using a well-established algorithm in the field of chaotic dynamics 33|. The 
idea is to evolve two initial conditions that start from a very close distance for about one 
dynamical time, then renormalize to bring the two particles close together again, and repeat 
the process until the average exponent associated with the orbital separation converges to 
an almost stable value. Typically convergence was achieved within ~ 100 orbital periods. 

After computing the orbits, we extracted the power spectrum for each orbit using a fast- 
Fourier-transform algorithm 32j. In doing so, we recorded each orbit at a rate ~ 40 times 
per orbital period. From the spectrum we computed the total power. Then we sorted the 
spectral frequencies in descending order, and starting from the highest frequency we added 
as many frequencies as were needed to reach 90% of the total power. The required number 
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of frequencies is defined to be the "complexity" n of the orbit 



3- 



2. Criterion for chaos 

Our first and foremost interest is to determine how many of the 2000 orbits in our sample 
are chaotic in a given TE configuration. Accordingly an objective, quantitative criterion 
for chaos is needed. There is no universally accepted criterion; hence, we developed our 
own using the following rationale. Both the largest short-time Lyapunov exponent \ and 
the complexity n are well-established, conventional measures of chaos 3^l- A first piece of 
information for defining the criterion comes from plotting n versus \ f° r a U °f the orbits. 
Fig. |21 provides an example. It shows that (a) n increases approximately linearly with Xi 
as one might expect since both quantities are measures of chaos, and (b) the regular orbits 
occupy a sharply defined region close to the origin of the n-vs.-x plot. The borders of this 
"region of regularity" thereby offer three possible criteria for chaos: one involving only the 
Lyapunov exponent, one involving only the complexity, and one involving both quantities. 

We chose to base our criterion only on the complexity for the following reason. Though 

no problem arises in computing Lyapunov exponents in the zeroth approximation, such is 

not the case in the exact potential, wherein we found that longer integration times are 

required to achieve adequate convergence. Recall that the exact potential must be specified 

over a three-dimensional grid. Accordingly, interpolation errors and discontinuities between 

the cells can affect computations of Lyapunov exponents because they involve the distance 

between two initially nearby orbits, which is a local property that is sensitive to the grid size 

and the order of interpolation. By contrast, a computation of complexity, in that it involves 

the Fourier spectrum of an individual orbit, avoids reference to nearby orbits and is thus a 

global property influenced little by the grid size, a notion that we corroborated during the 

course of our numerical studies. For simulations in exact potentials, we chose a grid size 

and interpolation algorithm (cf. Sec. IVII below) such that numerical errors had negligible 

effect on the computation of individual orbits. A standard measure of the "goodness" of an 

I I 

orbital integration is the degree to which total energy is conserved [36]; for every orbit we 
achieved conservation of total energy with relative error < 10 -6 . This is some two orders of 
magnitude better than contemporary standard practice. 

Investigations of n-vs.-x plots for many zeroth-order TE potentials led us to a specific 
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quantitative criterion, namely, an orbit is categorized as chaotic if its complexity n > 20 [cf. 
Fig. 121(c)]. This choice was also checked and confirmed by carefully inspecting hundreds of 
plots of individual orbits in several geometries. 

Plots of two representative orbits, one regular and one chaotic, as well as their power 
spectra in the x-direction and their surfaces of section in the dz/dt-vs.-z phase space, appear 
in Fig. El Both orbits have similar total energies and evolve in $o corresponding to Case 
5 and a 2 = 1.0, c 2 = 0.25, a strongly oblate spheroid. The first orbit has n = 8 and 
X = 0.028 tp 1 ; the second orbit has n = 178 and x = 0.640 t^ 1 . Their differences are 
striking. One major difference concerns their power spectra. The regular orbit exhibits two 
very distinct frequencies in its spectrum, whereas the chaotic orbit features a near-continuum 
spanning a large number of frequencies. The surfaces of section are computed by recording 
dz/dt and z when x = y = 0. The regular orbit is seen to stay within a localized region 
of the dz/dt-vs.-z phase space, whereas the chaotic orbit largely fills a global area of phase 
space commensurate to its total energy. 

Results of a numerical experiment that highlights the largest Lyapunov exponent, i.e., 
the rate of global chaotic mixing, appears in Fig. |3J Here, four initially localized clumps 
comprising 500 chaotic orbits evolve in the zeroth-order potential $ corresponding to Case 
5 and a 2 = 0.5, c 2 = 1.5, a triaxial configuration. Each clump initially occupies a cube of 
size 0.05 3 in the configuration space. The figure reveals that each clump mixes through a 
global region of phase space with an e- folding time comparable to a dynamical time tjj, taken 
here to be the orbital period corresponding to the total energy of the individual particles 
comprising the clump. After some tens of to each clump has spread through a volume 
commensurate to the total particle energy. 

C. Survey results 

Our strategy for surveying the parameter space of the TE configurations is as follows. We 
conduct the survey using the zeroth-order potential $ found from Eq. (JT7J). Recall that this 
potential depends on the focusing strength fl and, through i?(x), the scale lengths (a, c); 
we thus choose a specific parameter set [a, c;Q]. Then we integrate 2000 initial conditions 
generated by uniformly sampling (1) the volume spanned by the configuration, and afterward 
(2) the volume spanned only by the Debye tail. Having established the criterion for chaos 
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(complexity n > 20), we identify and count the chaotic orbits in the respective sample, and 
we express this number as a percentage of the sample. This percentage may be viewed as 
an indication of the extent to which a given parameter set supports globally chaotic orbits. 
However, it should not be taken too literally in that the initial conditions are distributed 
uniformly through a volume; they are not weighted by the actual density distribution. 

Our numerical experiments fall largely into two categories. Category I pertains to keeping 
Q fixed to its Case 5 value (intermediate space charge), i.e., i = 5 in Eq. (|18)L and then 
varying a and c within the constraint of Eq. (|13|). Category II pertains to keeping c 2 = 0.5 
fixed, and then varying Q and a. 

Figures and El depict the percentage of chaotic orbits for a portion of the Category I 
experiments. Also shown in the bottom panels of these figures are the initial conditions, 
projected onto the (x, z)-plane, corresponding to the axisymmetric configurations for which 
a 2 = 0.5. In addition, results from analogous experiments in exact potentials are plotted for 
comparison. The figures exhibit a number of features. First and foremost is the indication 
that most of the configurations support a considerable population of chaotic orbits. This 
is true even for axisymmetric configurations, but of course it is not true for spherically 
symmetric configurations in that their potentials are integrable. Second, orbits for which 
the initial conditions all lie within the Debye tail reflect a higher percentage of chaos than 
were they distributed through the entire configuration space. Third, prolate axisymmetric 
configurations (for which only one datum is shown in each of these figures) support little 
chaos. Fourth, for reasons to be elaborated in Sec. EH below, the exact potentials support 
less chaos than their zeroth-order counterparts. 

Further analysis reveals that, for configurations in which they are present, essentially 
all of the chaotic orbits originate in the Debye tail. Figure |7] dramatically illustrates this 
finding. To assemble this figure, the initial conditions are sorted and binned into increments 
of R spanning 0.1 units of length. Then, in each increment, the complexity n of every orbit 
is computed, and the complexities are averaged to obtain (n). The process is repeated for 
all of the increments, and then for many different parameter sets. The results, (n(R)) versus 
R for Case 5 with a 2 = 0.5 and several choices of c 2 , comprise Fig. 

Figure |H1 depicts the percentage of orbits for the Category II experiments for which the 
initial conditions are uniformly distributed over the volume spanned by the configuration. 
Cases corresponding to intermediate space charge would seem to support more chaos. This 
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finding makes sense when juxtaposed against the limiting cases of zero space charge at 
one extreme and zero Debye tail, i.e., the uniform beam, at the other extreme. With zero 
space charge, only the linear forces of the external potential influence the particles. For 
the uniform beam the external and space-charge potentials cancel one another so that the 
particles move freely, apart from reflections at the boundary of the configuration. In both 
extremes the orbits are all regular, excepting billiard effects, if any, associated with shapes of 
boundary surfaces of uniform bunches. The figure also suggests that prolate, axisymmetric 
TE configurations support little chaos. 



VI. NUMERICAL EXPERIMENTS IN EXACT POTENTIALS 

As a matter of principle, one must be concerned about the extent to which subtle structure 
in the potential can influence the qualitative behavior, and in particular the chaoticity, of 
an orbit. One well-known example is that of the Toda potential; the full Toda potential 



a truncated Toda potential is 
] . Our survey of the parameter 



is integrable and supports only regular orbits, but general 
not integrable and supports a population of chaotic orbits 
space of TE configurations centered on the use of $o ; a generally crude approximation to 
the true potential. The survey suggests a large region of the parameter space supports 
sizeable populations of chaotic orbits, wherein all of these orbits reach into the Debye tail. 
We may expect in general that the density profile, particularly that of the Debye tail, 
corresponding to the exact potential is considerably different from that corresponding to 
$q. For example, in the limit of distances very far from the centroid, the exact space-charge 
potential will approach spherical symmetry, whereas $o is everywhere homeoidally striated. 
Accordingly, to check and have confidence in the qualitative results of Sec.|VJ we must repeat 
the numerical experiments in a suitably broad collection of exact potentials. As mentioned 
earlier, the reason we did not base the survey on exact potentials is that the respective 
numerical experiments are computationally expensive. 

To integrate Eq. (JT2j) governing the exact potential $ (x) , which is a fully three- 
dimensional partial differential equation (PDE), we chose a multigrid algorithm |37[. The 
algorithm requires boundary conditions be specified over the surface of the volume occupied 
by the grid. We chose a cubic grid volume greatly exceeding the volume of interest, i.e., 
that spanning the Debye fall-off of the density. Then we calculated the boundary conditions 
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over the surface of this volume using the formalism of the first approximation, specifically, 
Eq. (|17|). Because the resulting boundary conditions are only first approximations to the 
true boundary conditions, we checked our numerical solutions by varying the positions of 
the bounding surfaces of the grid by factors of five, and we found negligible change in the 
results over the volume of interest. Applying the multigrid algorithm to three dimensions 
involves nontrivial manipulations of the inherent restriction, interpolation, and relaxation 
routines. In the process, a nonlinear algebraic equation emerges due to the nonlinearity 
of the PDE. It was solved using an iterative method that combines Newton-Raphson and 
bisection techniques 

After tabulating the exact potential in three dimensions, we used a fifth-order Runge- 
Kutta algorithm with variable time step to evolve the individual orbits in time, within 
which the force at every time (thus position) step was computed using a three-dimensional 
interpolation scheme. The resulting orbits conserved total energy with relative error better 
than 10~ 6 and sometimes as low as 1CT 7 . 

FigurelHlexemplifies the difference between the zeroth-order and exact potentials. The top 
two panels present isopotential contours in the (x, z)-plane for Case 5 with a 2 = 0.5, c 2 = 1.5, 
a triaxial configuration. The bottom panels show how the two potentials compare along 
each of the (x, y, z)-axes. Obviously there are, and there should be, differences, but the 
important question is to what extent these differences alter the qualitative evolution of the 
orbits and, in turn, the complexities that characterize them? Analogous graphs for Case 
5 and a 2 = 1.0, c 2 = 0.25, a strongly oblate spheroid, appear in Fig. from which the 
respective differences are seen to be much more pronounced. 

Figure ^2 provides a visual comparison of a chaotic orbit starting from the same initial 
condition and evolving in the zeroth-order and exact potentials of Fig. El Although orbits in 
the two potentials differ quantitatively, in many cases they are qualitatively similar in that 
they explore a similar volume of phase space and have similar morphology. This pertains 
to the example of Fig. ^2 however, this one example does not in any way guarantee every 
orbit that is chaotic in the potential of the zeroth approximation is also chaotic in the exact 
potential. Statistical comparisons of orbital complexities respective to the zeroth-order and 
exact potentials for several Case 5 configurations appear in Fig. ^1 and these show that 
the complexities in the two potentials can differ considerably depending on the specific 
parameter set under study. 
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Following the procedure delineated in Sec. IV CI we also computed the percentage of 
chaotic orbits in a broad range of exact Case 5 potentials. The results, juxtaposed against 
their counterparts computed using $ > appear in Table I and in Figs. Eland El For these 
examples, there is generally a smaller percentage of chaotic orbits in the exact potential than 
in the zeroth approximation. The explanation is simple: compared to the density ni[i2(x)], 
the density derived from the exact space-charge potential $ (x) is quasi- uniform over a larger 
volume and falls to small values over a shorter scale length. Accordingly, the configuration- 
space volume over which the space-charge force is markedly nonlinear, i.e., the Debye tail, 
is smaller. Figure El illustrates the difference in the density profiles corresponding to the 
approximate and exact solutions. In the limit of spherical symmetry the profiles are identical, 
and they disagree more strongly as they become less spherically symmetric. Most notable 
is the comparison between Figs. fTBT e) and (f ) concerning a strongly oblate spheroid, where 
we see that the corresponding exact density distribution is much more uniform than that of 
the zeroth approximation. This accounts for the strong discrepancy revealed in Fig. fT2T d) 
concerning orbital chaoticity. It is also consistent with expectations based on first principles: 
the closer a system is to being one-dimensional (e.g., sphere, cylindrically symmetric disc, 
infinite symmetric cylinder, in which particle motion is integrable), the less is the population 
of chaotic orbits. However, the essential observation is that, in most cases, the exact TE 
configurations do indeed support substantial populations of chaotic orbits in keeping with 
expectations that surfaced from the survey based on approximate solutions. 

VII. SUMMARY, IMPLICATIONS, AND FUTURE WORK 

We have explored orbital dynamics and phase mixing in thermal-equilibrium beams for 
which the potential is the superposition of an external potential quadratic in the coor- 
dinates and the self potential arising from space charge. The associated parameter space 
spans the full range of symmetries, i.e., spherical, cylindrical, and triaxial, and the full range 
of density profiles, i.e., gaussian (corresponding to negligible space charge) through uniform 
(corresponding to maximal space charge). To reiterate, the main findings concerning chaos 
in these systems, "discovered" in the context of zeroth approximations to the space-charge 
potentials and affirmed with the respective exact potentials, are: (1) configurations cor- 
responding to a large portion of the parameter space support considerable populations of 
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chaotic orbits, (2) essentially all of the orbits that are chaotic reach into the Debye tail where 
the collective space-charge force is manifestly nonlinear, (3) prolate axisymmetric configura- 
tions support little chaos, but prolate triaxial configurations can support considerable chaos, 
and (4) strongly oblate spheroids support little chaos, but moderately oblate spheroids can 
support considerable chaos. 

It is of interest to compare theoretical predictions concerning TE configurations, for which 
we herein have established the existence of chaotic orbits, with results of our numerical 
experiments. In terms of the dimensionless quantities introduced in Sec. IIII[ the parameters 
K and p of Eq. © take the form 



1 

K — - 

2 



1 1 



^(-j + l + ? )-(-(x)> 



yV(x) - (n(x)>a> 



P=1 ^ • (19) 



and Eq. (jHJ) then yields the mixing rate x- A comparison between theory and numerical 
experiments appears in Fig. El wherein the simulation results reflect statistics from initially 
localized clumps of 2000 particles that were started at zero velocity at various points in 
configuration space corresponding to various total particle energies E. The figure presents a 
plot of the mixing rate \ versus \E\ in the Case 5 configuration with a 2 = 4/5, c 2 = 4/3, a 
slightly triaxial system. This configuration is "not too far away" from spherical symmetry, 
which means the zeroth approximation $(x) = $o is correspondingly reasonable. It also 
means only a modest population (~ 5% for this parameter set) of chaotic orbits is supported. 
The figure was derived within the framework of the zeroth approximation because therein 
the Lyapunov exponents, i.e., mixing rates, can be accurately computed from the simulations 
and the microcanonical averages required for the theory likewise can be easily and accurately 
evaluated. The numerical experiments span a range 0.5 < \E\ < 60, corresponding to 
11 < R < 25, i.e., extending from within to well beyond the Debye drop-off in the density 
profile. The agreement between theory and numerical experiments is remarkably close. 

One can see from the numerical experiments described herein that chaotic mixing takes 
place on an e- folding time scale comparable to a dynamical time (an orbital period). This 
is very fast compared to, e.g., collisional relaxation; hence, one must account for this col- 
lisionless process when designing an accelerator for the production of high-peak- current, 
high-brightness beams. For example, particles comprising a beam out of equilibrium will, if 
globally chaotic, redistribute themselves globally and irreversibly on a dynamical time scale. 
Because perturbations induced, e.g., by transitions in the beamline will drive a beam away 
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from equilibrium, chaotic mixing can be a dynamic of practical importance. Consider the 
case of a TE configuration: a small perturbation from image charges passing through an 
external irregularity in the beamline will distort the Debye tail. If a substantial fraction of 
particles in the Debye tail are chaotic, which is the case for a wide range of bunch geometries, 
a corresponding fraction of the orbits comprising the distortion will quickly mix throughout 
the volume of the configuration. The work done by the external perturbation in setting up 
the distortion will thereby appear in the form of a larger configuration-space volume. If the 
perturbation is strong enough so that mixing in momentum space associated with conse- 
quent time-dependence in the potential is also substantial, then some of the work done will 
also appear in the form of a larger momentum space. The net effect is a larger emittance. 
If there are many such perturbations along the beamline, the cumulative emittance growth 
may be troublesome. 

The present investigation and its associated implications concern only very specific, time- 
independent, single-species systems, i.e., beams (or nonneutral plasmas) in thermal equilib- 
rium. These are the most benign systems imaginable, yet we found even they can support 
chaotic orbits. Any perturbation will create a nonequilibrium, time-dependent system that 
will subsequently evolve self-consistently. Accordingly, the space-charge potential can be 
complicated, particularly if the perturbation is strong. The only sensible conjecture un- 
der such conditions is that the corresponding population of chaotic orbits will be larger, 
and in turn chaotic mixing will be more prevalent. Exploratory numerical simulations of 
an equipartitioning system and of merging beamlets have supported this notion Fur- 
ther exploration of time-dependent beams is warranted and will likely prove illuminating, 
particularly in regard to deciphering time scales for emittance growth, halo formation, etc. 

By using only smooth potentials we have restricted our analysis to the six-dimensional 
phase space of a single particle. Accordingly we have suppressed dissipative effects of colli- 
sions in particular, and force fluctuations in general. Such effects can only enhance chaos, 
as has been demonstrated, e.g., in numerical experiments concerning self-gravitating sys- 
tems I22I. As the next step, we have constructed frozen iV-body representations of the 
charge densities of the TE configurations and with these representations are repeating the 
numerical experiments described herein. One of our objectives is to determine the mini- 
mum number of particles needed to reproduce the dynamics associated with smooth time- 
independent potentials. Results will be described in a forthcoming paper [28 1. In the future 
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it will be of interest to do likewise for time-dependent systems and ultimately ascertain, e.g., 
conditions under which the Vlasov equation governing the six-dimensional phase space of a 
single particle can be applied with confidence. 
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Table I. Percentage of chaotic orbits: approximate vs. exact potentials 



Whole configuration space: 



a 1 


c 2 


% th 


% exact 


0.5 


0.50 


6.95 


5.90 


0.5 


1.50 


35.50 


19.45 


0.5 


2.00 


36.55 


17.15 


1.0 


0.25 


33.20 


7.05 


1.0 


0.50 


19.85 


26.05 


Debye fall-off interval: 


a 2 


2 

C 


% th 


% exact 


0.5 


0.50 


7.55 


7.32 


0.5 


1.50 


43.40 


35.50 


0.5 


2.00 


42.90 


36.55 


1.0 


0.25 


37.55 


7.05 


1.0 


0.50 


29.45 


24.55 
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FIG. 2: (a) Distribution f(x), in arbitrary units, of Lyapunov exponents x (Case 5, a 2 = 0.5, c 2 = 
1.5) in the zeroth approximation $o- The unit of x is to'- The distribution peaks at low values of 
X for which the respective orbits are regular. The dotted line, hand-drawn at the point where the 
distribution levels off, suggests one possible criterion of chaos, (b) Distribution f(n) of complexities 
n corresponding to the orbits of the first panel, which looks qualitatively similar to /(x)- (c) 
Complexities n versus Lyapunov exponents x- The inset reveals that the concentration of regular 
orbits near the origin lies inside sharply defined boundaries. 
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FIG. 3: Two representative orbits in the potential $o of the zeroth approximation (Case 5, 
a 2 = 1.0, c 2 = 0.5), one regular (a,c,e) for which n = 8 and x = 0.028 , and one chaotic (b,d,f) 
for which n = 178 and x = 0.640 tp, each with similar total energies: (a,b) orbits on (x,z)-plane; 
(c,d) power spectra in x-direction; (e,f) surfaces of section dz/dt vs. z. 
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FIG. 4: Snapshots of four different clumps of chaotic orbits (Case 5, a 2 = 0.5, c 2 = 1.5) evolving 
in the potential $o- Each clump initially occupies a cube of volume 0.05 3 , but exponentially grows 
to fill a volume commensurate to the total particle energy. 
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c 2 = 0.25 c 2 =0.50 c 2 =1.00 c 2 =1.50 c 2 =2.00 



FIG. 5: Top: Percentage of chaotic orbits for Case 5 vs. c 2 in the potential $o with different 
choices of a 2 . The initial conditions uniformly sample the region < R < 15, which in essence 
covers the volume spanned by Case 5 configurations [cf. Fig. IHb)]. Results derived from exact 
solutions ^(x) of Poisson's equation are also plotted for comparison; they are joined by a dotted 
line (a 2 = 1, c 2 < 1) or by dashed lines (a 2 = 0.5). These exact results are to be compared to 
their zeroth-order counterparts delineated with blackened symbols to aid the eye. Bottom: Initial 
conditions, projected onto the (x, z)-plane, used for a 2 = 0.5; z-axis is vertical. 
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FIG. 6: Same as Fig. [51 but with initial conditions that uniformly sample only the Debye tail 
9 < R < 15. For ease of visualization, the bottom panel shows only the initial conditions for which 
V ~ 0. 
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FIG. 7: Average complexity (n) versus homeoidal coordinate R for Case 5 in the potential <l>o 
with a 2 = 0.5 and different choices of c 2 . Essentially all of the chaotic orbits reach into the Debye 
tail. 



32 




FIG. 8: Percentage of chaotic orbits in the potential <l?o with c 2 = 0.5 and Q corresponding to 
Cases 3-8 as defined in Eq. (|18|), Different curves correspond to different choices of a 2 . The initial 
conditions uniformly sample the respective configurations. 
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FIG. 9: Comparisons between the space-charge potential <£o m the zeroth approximation and the 
exact potential 3>(x) (Case 5, a 2 = 0.5, c 2 = 1.5): (a) isopotential contours of (b) isopotential 
contours of 3>(x), (c,d,e) profiles along the (x, y, z)-axes, respectively (blue and red curves pertain 
to $0 and 3>(x), respectively). 
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FIG. 11: A chaotic orbit in the zeroth-order potential (top panels) and exact potential (bottom 
panels) of Fig. ED evolved from the same initial conditions. The orbit is similar, but not identical, 
in the two potentials. 
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FIG. 12: Average complexity (n) versus homeoidal coordinate R for zeroth-order (blue) and exact 
(red) Case 5 potentials with: (a) a 2 = 0.5, c 2 = 0.5; (b) a 2 = 0.5, c 2 = 1.5; (c) a 2 = 0.5, c 2 = 2.0; 
(d) a 2 = 1.0, c 2 = 0.25; (e) a 2 = 1.0, c 2 = 0.5. 
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FIG. 13: Comparisons between the density distributions of the first approximation ni[i2(x)] (left) 
and exact solution n(x) (right) for Case 5 with: (a,b) a = c = 1.0 (spherical), (c,d) a 2 = 0.5, c 2 = 
1.5 (triaxial), (e,f) o 2 = 1.0, c 2 = 0.25 (strongly oblate spheroid). 
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FIG. 14: Theoretical results (dashed curve) and numerical results (diamonds) for the mixing rate 
X of chaotic orbits vs. total particle energy \E\ (Case 5, a 2 = 4/5, c 2 = 4/3). The unit of x is '• 
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